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Abstract 

Recently, a novel strategy for post -experiment state estimation of discretely-measured dy- 
namic systems has been developed. The method accounts for errors in the system dynamic model 
equations in a more general and rigorous manner than do filter-smoother algorithms. The dynamic 
model error terms do not require the usual process noise assumptions of zero-mean, symmetrically 
distributed random disturbances. Instead, the model error terms require no prior assumptions 
other than piecewise continuity. The resulting state estimates are more accurate than filters for 
applications in which the dynamic model error clearly violates the typical process noise assump- 
tions, and the available measurements are sparse and/or noisy. Estimates of the dynamic model 
error, in addition to the states, are obtained as part of the solution of a two-point boundary value 
problem, and may be exploited for numerous reasons. In this paper, the basic technique is ex- 
plained, and several example applications are given. Included among the examples are both state 
estimation and exploitation of the model error estimates. 
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1.0 Introduction 


A large number of applications exist in the general area of “post- experiment” or “post- 
flight” estimation, wherein estimates of the state histories of a dynamic system are obtained using 
an assumed state dynamic model and sets of discrete measurements. In general, both the assumed 
model and the available measurements are imperfect. The motivation for applying an “optimal 
estimation” algorithm is to combine the model output with the available measurements in such 
a way as to obtain estimates of the state histories which are superior to both the model and the 
measurements, and, in addition, satisfy an optimality criterion. In this paper, a new estimation 
strategy is described which includes both a new optimality criterion and a new algorithm for 
obtaining estimates based on this condition. 

The following generic problem statement for post-experiment estimation of a dynamic pro- 
cess is used to motivate the discussion. Given a system whose state vector dynamics is modeled 
by the (linear or nonlinear) system of equations, 

* = ./[*(*)»«(*)»*] (i) 

where 

x = n x 1 state vector 
f = n X 1 vector of model equations 
u = p X 1 vector of forcing terms , 

and given a set of discrete measurements modeled by the (linear or nonlinear) system of equations, 


£( tfc ) = £*[£(**)> 4 k] + v k , k = 1 (2) 

where 

y(tk) ErXl measurement set at t*. 

^Erxl measurement model equations 
m = total number of measurement sets 
y_ k = r X 1 measurement error vector , 

and Vk is assumed to be a zero-mean, gaussian random sequence of known covariance Rk, determine 
the optimal estimate for x(t) (denoted by i(t)), during some specified time interval t Q < t < tf . 
Clearly, the definition of optimal is subjective, and we begin by discussing optimality criteria. 

2.0 Optimality Criteria 

The typical approach for obtaining an optimal estimate of the system state trajectories is 
the minimization of a function of the estimate error, 

£* = £{(4 -a)} (3) 


or its covariance, 

= (4) 

Among these criteria are the well-known “maximum likelihood” and “minimum variance” strategies 
(e.g., Gelb 1 ). For example, the minimum variance criterion requires the minimization of the trace 
of j P ii. Many other criteria which rely on estimating the estimate error statistics (Eqs. (3), (4)) 
have also been used as bases for estimation algorithms. A practical problem arises during actual 
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implementation of these methods. In order to determine the estimate error statistics, it is necessary 
to assume that any errors in the system model Eq. (1) are noise of known probability distribution. 
Most often, the probability distribution is assumed to be zero-mean gaussian, whose covariance 
is treated as a known quantity (“process noise”). The state estimation proceeds without any 
adjustment to the system model equations. 

In general, it is difficult, if not impossible, to rigorously justify process noise assumptions for 
the model error. For real physical systems, error is likely to be due to modeling simplifications such 
as linearization, neglect of higher-order terms, etc., or, perhaps, just plain ignorance. Many of these 
likely sources are deterministic and non-zero mean. Consequently, the estimate error statistics in 
Eqs. (3) and (4) cannot be calculated rigorously. Estimates based on their optimization are sub- 
optimal, e.g., the minimum variance estimate is not truly minimum variance if the variance which 
is minimized is not the true variance. 

These observations are well-known and are repeated here only to motivate the discussion. 
The optimal estimation strategies which require process noise assumptions work well in many 
applications, whether or not the model error assumption is justifiable, and filter algorithms are 
the most commonly used estimators in practice. The filters must generally be artistically “tuned”, 
but this is often possible and sufficient for a reasonably accurate estimate. 

However, filter accuracy may deteriorate substantially under a number of conditions. The 
filter algorithms rely on the integration of the original dynamic model Eq. (1) for the between- 
measurement estimate. If the model is poor and the measurements are sparse, the accumulated 
integration error between measurements may become very large. Even if the measurements are 
dense, if they axe particularly noisy, and the model is poor, then the filter estimate may be of 
poor accuracy. Under certain conditions, the filter may become unstable. Divergence of filters 
when process noise assumptions are violated may be found in Fitzgerald 2 , Huber 3 , and Breza and 
Bryson 4 , among others. 

With this motivation, Mook and Junkins 5 developed a new estimation strategy which elim- 
inates any apriori assumptions about the model error except that it is continuous between the 
measurement times. The method, called Minimum Model Error (MME) estimation, is based on an 
optimality criterion which does not require estimation of the estimate error statistics, Eqs. (3)-(4). 
In the remainder of this paper, a summary of the method is given, followed by several application 
examples. 


3.0 The Covariance Constraint Optimality Criterion 

In the MME, a novel optimality criterion is used. The probability distribution of the state 
estimate error is not estimated. Instead, the optimal state trajectory estimate is determined on 
the basis of the assumption that the measurement-minus-estimate error covariance matrix must 
match the measurement -minus- truth error covariance matrix. This condition is referred to as the 
“covariance constraint”. The covariance constraint is defined mathematically by requiring the 
following approximation to be satisfied: 

tj) ][y(^j) ~ 9_ | 555 Rj (5) 

Thus, the estimated measurements g(i{tj) y tj) are required to fit the actual measurements y(tj) 
with approximately the same error covariance as the actual measurements fit the truth. An algo- 
rithm for obtaining the estimates is described shortly. 


133 


The covariance constraint may be evaluated without knowledge of the estimate error statis- 
tics. Consequently, there is no need for process noise-like assumptions for the model error. In 
the next section, an algorithm which produces estimates which satisfy the covariance constraint is 
derived, treating model error as an unknown which is estimated along with the states. 

The interpretation of the “approximately equal” sign in the covariance constraint may be 
adjusted according to the particular application. If the measurements are repeated samples of 
the same quantities, as is usual in a filtering problem, then a good approach is to calculate the 
covariance of the measurement-minus-estimate residuals using all of the measurements simultane- 
ously. Thus, the covariance constraint is averaged over all of the measurements. In problems where 
several distinct sets of measurements are repeated, each set may be averaged separately. An exam- 
ple is spacecraft navigation, where the measurements may include sets of attitude measurements 
and sets of angular velocity measurements. These two sets are normally taken independently and 
contain different noise levels, so they should be averaged separately. 

4,0 MME Algorithm 

If the dynamic model Eq. (1) contains significant error, then its output generally cannot 
predict the measurements with enough accuracy to satisfy the covariance constraint. The estimated 
measurement set at time tk is based on the current state estimate, £(£*), as shown in Eq. (2). The 
between- measurement state estimate is based on integration of the system dynamic model. Thus, if 
the system dynamic model contains errors, the integration does not yield the correct state estimate, 
and the residuals between the estimated and the actual measurements are too large. Consequently, 
the model error must be reduced in order to satisfy the covariance constraint. To accomplish this, 
a model correction term d(t) is added to the original dynamic model as 

*= j&(t),u{t),t) + d{t) (6) 

In general, an infinite number of d(£)’s exist which are capable of correcting the model to satisy the 
covariance constraint. The minimum correction is sought, thereby providing the least adjustment 
to the original model. Accordingly, the following cost functional is minimized with respect to d(t): 

i 1 

f tf rj, 

+ / d T {r)Wd(r)dr (7) 

Jt 0 

where W is a k X k weight matrix chosen to satisfy the covariance constraint as described shortly. 
The functional J in Eq. (7) is the sum of two penalty terms, whose relative weighting is controlled 
by W . If W is near zero, then the integral term is nearly zero. Consequently, the allowable 
d(t) is virtually uhUmited and thus the model is corrected until the measurements are predicted 
almost exactly (i.e,, the summation term goes to zero). However, this is only appropriate when the 
measurements are perfect. If the measurements are noisy, then the covariance constraint implies 
that the summation term should not be zero. The weight matrix, W , is chosen such that the 
covariance constraint is satisfied, allowing just enough correction d(t). Generally, determination of 
W requires a search procedure. However, unlike the tuning of a filter, which is essentially artistic, 
W is specified by satisfaction of the covariance constraint. 

In Figure (1), the concept of the covariance constraint is demonstrated for a one- dimensional 
(n=l) system. Figure (la) shows the ratio of the left-hand side of Eq. (5) to the right-hand 
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Figure 1. Choosing W to satisfy the covariance constraint leads to the optimal state 
estimate. 

side, plotted versus W. As the weight is decreased, the corrected model predicts the actual 
measurements more closely as shown in the figure. In Figure (lb), the estimate variance is plotted 
versus W . The covariance constraint implies that the optimal estimate occurs when the covariance 
constraint is satisfied. 

An algorithm for the minimization of J in Eq. (7) follows directly from a modification (e.g., 
Geering 6 ) of the so-called Pontryagin’s necessary conditions (e.g., Rozonoer'). For a given IF, the 
minimization of J in Eq. (7) with respect to d{t) leads to the two-point boundary value problem 
(TPBVP) summarized as: 


i ~ f[x(t),u(t),t] + d(t) 
\dx/~ 


d — ~-W~ l 

2 


-di 

du 


iT 


( 6 ) 

( 8 ) 

( 9 ) 
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x(f-o) = specified , or A (t 0 ) = 0 

(10) 

M*j ) = M f J ) + 2 Hf RJ 1 [y(tj ) - £(*(*;),*;)] 

(11) 

x(t j ) = specified , or A(t+ ) = 0 

(12) 


where 



)»*> 


This TPBVP contains jump discontinuities in the costates at each measurement time where the 
predicted measurement does not exactly match the actual measurement. The size of the jump 
is proportional to the measurement residual \y{tj) — £(i(tj),tj)], which, via the covariance con- 
straint, is proportional to the measurement noise. Prom Eq. (9), these costate jumps lead to 
jumps in the estimated model error. Thus, for noisy measurements, the model error estimates are 
jump discontinuous proportional to the measurement residuals. Note that this is identical to a 
filter except that in a filter, the jumps are in the state estimates. The MME state estimates are 
continuous. 

The algorithm Eqs. (6)-(12) exhibits several desirable features of both batch and sequential 
estimation techniques. The state estimate is obtained by processing all of the available measure- 
ments, much like a batch estimator such as least squares. Thus, the estimate is optimized in a global 
sense. In addition, the state estimate is continuous, eliminating the jump discontinuities present 
in filter estimates. For many physical systems, jump discontinuities in the states are not possible; 
thus, jump discontinuities in the filter state estimates must be reconciled in an artful manner. In 
addition to the batch algorithm-like advantages, the minumum model error algorithm calculations 
are based upon sequential processing of the measurements, which, like the filter algorithms, greatly 
reduces the memory requirements and eliminates the need for large matrix manipulations. From 
the standpoint of algorithmic calculations, the minimum model error technique shares advantages 
of both batch and sequential estimation techniques. 

If the assumed model in the MME algorithm is linear, then a multiple shooting technique 
may be used to solve the TPBVP described by Eqs. (6)-(12) (Lew and Mook 8 ). This technique 
converts the TPBVP into a set of linear algebraic equations which may be solved using any linear 
equation solver. 

When the covariance constraint has been satisfied, the estimate is considered to have been 
optimized. As a byproduct of the solution, the estimates d(t) of the model error required to 
satisfy the optimality criterion are available. The results of the examples clearly indicate that 
these terms may provide highly accurate estimates of the actual model errors, leading to potential 
improvements in the model. 


5.0 Examples 

In this section, several example applications are summarized which demonstrate the present 
method and explore the accuracy of both the state estimates and the model error estimates obtained 
using it. The examples include both linear and nonlinear systems, varying degrees of model error, 
varying levels of measurement noise, measurement frequency, and total number of measurements. 
Exploitation of the model error estimates is also demonstrated. 
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5.1 Simple Example of Minimum Model Error Estimation 

To illustrate the application of the minimum model error approach, consider estimation of a 
scalar function of time for which noisy measurements are the only information available. No prior 
knowledge of the underlying dynamics is assumed. Thus, the system dynamic model equation is 

x = 0 (13) 


Using the minimum model error approach, the system model is modified by the addition of a 
to-be-determined unmodeled effect as 

x = 0 + d(t) (14) 

where d(t) represents the dynamic model error. For simplicity, the measurements are direct mea- 
surements of the state itself, and the measurement noise is a zero mean gaussian process with a 
variance of a 2 , given as 

x(t k ) = x{t k ) + v k , k ~ 0 , m (15) 

where x k is the measurement at time tfc, Xk is the true state at time t k , and v k is a zero-mean 
gaussian sequence of variance a 2 . The cost functional to be minimized (see Eq. (7)) is 


i m r tm 

j = — - *(**)]* + / d2 ( r ) Wd 

a fc =0 Jt ° 


(16) 


where W is the to-be-determined weight on the integral sum-square model error term. The TPBVP 
which results from the minimization of J with respect to d(t) may be summarized as (see Eqs. (6)- 
( 12 )) 


d= - 


2 W 


x = d = — 


2 W 


A = -^A = ° 

OX 

Kt o") = A(t+) = ° 

fc ) + ~^2 ~ 

where A is the vector of costates. The boundary conditions indicate that the state is unknown at 
to or £/. The algorithm proceeds according to the following steps: 


1) Choose W 

2) Solve the TPBVP 

3) Check the covariance constraint 

4) If the covariance constraint is not satisfied, go to step 1 

The true state history for this example is taken as x(t) = cos(t). In Fig. 2, a set of 101 
simulated measurements spanning the time interval t 0 — 0 to tf = 10 is shown. The measurements 
were simulated by adding a computer-generated gaussian random sequence to the true state as 


x k = cos(t k ) + Vk 


(17) 
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Figure 2. Simulated measurements of cos(t) with <r 2 = 0.114. 



Figure 3. MME estimates using the Figure (2) measurements with no model. * denotes 
measurement, denotes truth (cos(t)), □ denotes MME estimate. 

The nominal variance of Vk in Fig. 2 is 0.1, although the actual variance depends on the seed 
supplied to the random number generator. This variance is 10% of the peak amplitude. Thus, the 
average measurement error is approximately 50% of the average amplitude. 

In Fig. 3, the minimum model error state estimate is shown along with the measurements 
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and the true state history. Note that the state has been reconstructed to an error variance of 
.00S5, considerably better than the measurement, variance even in the total absence of a model. 
Note also that the model prediction variance (i.e., constant x — 0) is 0.717. 



Figure 4. MME model error estimates using the Figure (2) measurements with no 
model. <> denotes true model error (— sm(t)), □ denotes MME estimate. 


In Fig. 4, the model error term is plotted along with the true model error, - sin(t ). Although 
the model error estimate contains considerable noise, due to the noisy measurements, it is an 
accurate representation of the actual model error. Based on an examination of Fig. 4, a user might 
easily conclude that the dynamic model error is indeed — If the dynamic model is amended 

from x = 0 to x = — sm(t), and the estimation process repeated, the state estimate is virtually 
exact and the model error estimate is virtually zero. 

5.2 System State Estimation from MME 

Several applications examples are now presented for system state estimation using the MME 
method. Other examples have also been investigated but are omitted here due to space limitations. 

5.2.1 Modal Space State Estimation 

In these examples, taken from Mook and Lin 9 and Lin 10 , the state histories of the output 
measurements of a system described by a linear sum of system modes are obtained using the MME. 
The simulated measurements are created by assuming a truth as a sum of several modes, and then 
adding gaussian noise to the truth. The assumed model for the MME estimation is taken as the 
first mode in the sum. All of the other modes are ignored in the assumed model. Figures (5) shows 
the results from a case with a five mode truth. The one-mode model is plotted along with the 
MME estimate and the truth. The model is seen to be very poor, but the estimate is essentially 
perfect. The measurement noise is gaussian with a 2 = 0.04, and the measurement interval is 0.1 
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Figure 5. One-mode assumed model, five-mode truth, and MME estimate (*). 
seconds. 

Clearly, the MME is able to recover from a very poor model to produce accurate state 
estimates. This result can be very helpful for structural modelers who axe uncertain about whether 
or not a modal model has been truncated with the correct number of retained modes. 


5.2.2 Nonlinear State Estimation 


The following example is taken from Mook 11 . Consider the single-degree-of-freedom system 
modeled by 

x + ulx = f{x(t),x{t)) + F(t) (18) 

where x and x are the system states, F(t) is a known external excitation, and /(z(t),£(£)) contains 
terms which may be nonlinear in the states. The external excitation is assumed to be independent 
of the states. Eq. (18) may be converted to state-space form as 


±(t) 


0 

m 


+ 


f(x{t),x(t) 


(19) 


where z = {x(t) x(t)} T . A specific example, after Thompson and Stewart 12 , is given by 


2.56.1; + 0.32.£ + x -f 0.05:c 3 = 2.5cos(t) 


( 20 ) 


This example exhibits two distinct possible steady-state solutions, depending on the initial con- 
ditions. The assumed model for the MME algorithm is used in two different forms. First, the 
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nonlinear system is modeled for the MME as a linear oscillator. Measurements are simulated with 
a variety of different noise levels and frequencies. In each case, the MME is able to obtain accurate 
state estimates. Results are shown in Figures (6) through (8). In part (a) of each figure, the 
measurements are plotted along with the linear model output, thus showing the information given 
to the MME method. In part, (b) of each figure, the truth, measurements, and MME estimate 
are shown. The accuracy of the state estimate is apparent from the figures, even for measurement 
frequency less than 4/cycle and total measurements as low as 15 (Figure (7)), and for noise levels 
with a equal to 30% of the peak amplitude (Figure (8)). 



Figure 6. (a) Nonlinear estimation with 30 noiseless measurements, using assumed 
linear oscillator model, (b) Truth, measurements, and MME estimates. 



Figure 7. (a) Nonlinear estimation with 15 noiseless measurements, using assumed 
linear oscillator model, (b) Truth, measurements, and MME estimates. 


Second, the assumed model for the MME consisted only of the external forcing, so that no 
knowledge of the system is assumed. Results are shown in Figure (9), where the MME estimate 
pictured in part (b) is seen to be very accurate despite the very poor model pictured in part (a). 
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Figure 8. (a) Nonlinear estimation with 100 noisy measurements, with noise level 
approximately 30%. (b) Truth, measurements, and MME estimates. 



Figure 9. (a) Nonlinear estimation with 100 noiseless measurements, using no assumed 
model, (b) Truth, measurements, and MME estimates. 



5.3 System Identification From Model Error Estimates 

An area of considerable interest in many engineering disciplines is identification, the process 
of obtaining an accurate model of a dynamic process using measured data. State estimation and 
identification are most often two separate processes. Some versions of Kalman filters have been 
implemented which treat unknown constant parameters in a model as states, so that they are 
estimated as part of the state vector. This approach, like most other identification techniques, 
requires the user to construct a model of appropriate form and order apriori. The filter then 
estimates the constant parameters in the model. However, the filters still assume that any model 
error is a gaussian white noise, so this approach usually works well only if the model order and 
form are correctly chosen by the user. The MME, by estimating the model error, may be used to 
determine the form of the model error before attempting to estimate the parameters. 

Several studies have been conducted to investigate the use of the MME as an aid to system 
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Figure 10. Armature-controlled motor drives a rotating shaft assembly with inertia 
J, damping 1 5, and stiffness k. The motor constants are Kt and Kb, and z = {0 0 } T . 


identification. These results are presented next. 

5.3.1 Linear State Space Parameter Identification 

The simplest form of model to identify is the linear, time-invariant state- space model. 
This is also the most commonly used model form in practice. Some examples of linear system 
identification axe given in Mook, Liu, and Ho 13 , and some results of that study are repeated 
here. Two assumed true systems are studied; an armature-controlled motor system driving a 
rotating shaft assembly, shown in Figure (10), and a linear, two-degree-of-freedom model nominally 
represented by two masses, springs, and dampers, as shown in Figure (11). In Figure (10), the term 
A 22 = ~ — 1 , which represents damping, is assumed to be unknown and is to be estimated. 

In Figure (11), the three damper constants are to be estimated from the free response. Simulated 
measurements are created for both cases by adding gaussian white noise to the truth, and several 
cases axe presented for varying noise, measurement frequency, and record length. In addition, the 
assumed model used for the MME is varied from case to case by altering the assumed values for 
the unknown constants. 

The parameter estimation is carried out by a least- squares fit of the estimated model error. 
Since the model error estimate is continuous except at the measurement times, it may be sampled 
at an arbitrary number of points away from the measurements to create an over determined system 
of algebraic equations in the unknown parameters. Then, a least- squares algorithm is used to 
produce the parameter estimates. Results for the system in Figure (10) are shown in Table (l), 
and for the system in Figure (11), in Table (2). 

5.3.2 Nonlinear System Identification 

In section 5.2.2, results are given which demonstrate very accurate state estimation of a 
nonlinear system, given poor dynamic models and noisy, sparse data. In this section, identification 
results are given for those same examples. More detail is available in Mook 11 . 

The model error estimates corresponding to Figures (6)-(9) are shown in Figures (12)- 


143 




*- Fi 



* F2 

♦ z2 
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-(*i+* 2 ) -(ci+C 2 ) -fc 2 -C 2 

0 0 0 1 
V -fc 2 -C 2 -(^2 + ^3) -(c 2 + c 3 )/ 


z + 



Figure 11. Mass-spring-damper system, where F is applied force. 


Table 1. Parameter estimates for the system in Figure (10). 


Measurement 
Frequency 
<{neas/10 sec) 

Measurement 

Variance 

True A Z 2 

Estimated A, 2 

Error 

41 

0* 

- 0. 4 

-0. 4000 

0. 00 

41 

0* 00018 

-0. 4 

-O. 394691 

1. 33 

21 

0. 00014- 

• 

0 

1 

-0. 395293 

1. 18 

11 

0. 00015 

-0. 4 

-0. 415468 

3. 87 

41 

0» 0045 

“0# 4 

-0.361596 

9.60 


(15). Clearly, the model error estimates are dependent on the accuracy and frequency of the 
measurements. For more accurate and more frequent measurements, the model error estimates are 
smoother and more accurate. However, the accuracy of the model error is not dependent on the 
accuracy of the assumed model. This is a very significant result for identification. 

In order to identify the nonlinear model from the model error estimates, a parameterized 
model of the model error is constructed and then the parameters are estimated using a least- 
squares algorithm. For demonstration purposes, the assumed model for identification contained 
more terms than the actual model error, including the case when no prior model is assumed for 
the MME (Figures (9) and (15)). The least squares fit produced near-zero parameter estimates for 
the assumed model error terms which are not in the model, and near-perfect parameter estimates 
for the assumed terms which are in the model. 
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Table 2. Parameter estimates for the system in Figure (11). 



(values are in the order of cl, c2, c3> 



Figure 12. Model error estimation Figure 13. Model error estimation 

with 30 noiseless measurements, with 15 noiseless measurements, 

using assumed linear model. using assumed linear model. 

5.2.3 Modal Space Realization/Identification 

Recently, considerable interest in the identification of modal space models for large flexible 
systems has arisen in conjunction with such proposed projects as the space station. Several methods 
which produce accurate modal models from time-domain data have been developed (e.g., Ibrahim 
and Mikulcik 14 ; Rajaram and Junkins 15 ; Hendricks et a/ 16 ; Chen et al 1T ). The recently developed 
Eigensystem Realization Algorithm (Juang and Pappa 18 ) is particularly attractive because it first 
determines the model order and then estimates the model parameters. This alleviates a very serious 
drawback of most methods, which require that the model order be known aprion . However, Juang 
and Pappa 19 found that for high noise levels in the measurements, the ERA could not determine 
the correct number of modes in the model, and the parameter estimates for the model were of low 
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Figure 14. Model error estimation 
with 100 noisy measurements, 
using assumed linear model. 



Figure 15. Model error estimation 
with 100 noiseless measurements, 
using no assumed model. 



Figure 16. Algorithm flowchart for modal identification, 
accuracy. 

In Mook and Lew 20 , the MME method is used in conjunction with the ERA to produce 
an algorithm which is significantly less sensitive to noise. The algorithm may be summarized as 
(i) apply ERA to the original measurements, (ii) use the ERA-produced model and the original 
measurements in the MME to produce state estimates of the measurements, (iii) sample the MME- 
produced state estimates to create simulated measurements of higher accuracy than the original 
measurements, and (iv) apply ERA to the simulated measurements in order to realize/identify the 
correct model. The steps (ii)-(iv) may be repeated, since the MME will produce more accurate 
state estimates if a more accurate model is used. Consequently, if the first pass through steps (ii)- 
(iv) produces more modes than step (i), a second pass through steps (ii)-(iv) may further improve 
the accuracy of the realization/identification. The entire procedure is represented by the flowchart 
in Figure (16). 

The combined algorithm has been investigated for identification of the modes of a clamped- 
clamped beam. The true model is given by 

y(t) = l.Osm(f) + 0.05sm(2.76t) + 0.001sm(5.4t) (21) 

Measurement data was created with several noise levels, including a = 0.001, a = 0.003, a — 0.01, 
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cr — 0.05, and a — 0.1. Note that, the highest noise level corresponds to approximately 10% of 
the measurement amplitude, while the lowest noise level is approximately 0-1% of the amplitude. 
Moreover, the highest noise level is twice as high as the second, and 100 times as high as the third, 
modal amplitudes, while the lowest noise level is equal to the smallest modal amplitude. 


Table 3* Singular values from the ERA algorithm. 


<7 = 0.1 

a = 0.05 

<7 = 0.01 

a = 0.003 

o = 0.001 

II 

O 

o 

29.484 

29.468 

29.458 

29.457 

29.457 

29.456 

20.038 

20.030 

20.028 

20.028 

20.028 

20.028 

2.166 

1.675 

1.331 

1.279 

1.264 

1.257 

1.638 

1.199 

0.904 

0.862 

0.851 

0.845 

1.196 

0.593 

0.118 

0.035 

0.023 

0.026 

1.181 

0.590 

0.113 

0.034 

0.019 

0.022 

1.134 

0.566 

0.110 

0.033 

0.012 

10" 13 

1.124 

0.561 

0.108 

0.031 

0.011 

10- 13 

1.083 

0.541 

0.106 

0.030 

0.011 

10-13 

1.031 

0.516 

0.102 

0.029 

0.010 

10-' 3 

0.965 

0.483 

0.097 

0.029 

0.010 

10-13 

0.942 

0.471 

0.094 

0.028 

0.009 

10-13 


The model order is determined from the singular values of the singular value decomposition 
of H( 0), where H is the so-called “Hankel matrix’-. The model order is determined by the number 
of pairs of singular values between which there is a significant drop in magnitude. Table (3) gives 
the singular values in order from largest to smallest for each of the five noise levels. It appears that 
the model order is one, perhaps two, for a = 0.1, since the singular value pairs beginning with the 
second pair are approximately the same magnitude. Consequently, for noise levels of a = 0.1, the 
ERA method indicates one mode from the measurements. This seems intuitively reasonable since 
the level of noise exceeds the amplitude of modes 2 and 3. 


Table 4. Parameter identification from the ERA algorithm. 



Mode 1 

Mode 2 



Mode 3 

Noise 

Frequency 

Damping 

Frequency 

Damping 

Frequency 

Damping 

(7 = 0.1 


- 2.25 xlO " 3 

2.473 

- 2.410 

11.48 

WEBmm 

<7 = 0.05 

1.013 

- 1 . 37 xl 0" 3 

2.660 

- 0.967 

12.87 


cr =0.01 

1.003 

- 3 . 23 x 10 “* 

2.739 

- 0.048 

26.26 


o = 0.003 

1.001 

- 8.37 xlO - 4 

2.753 

2 . 11 x 10 -* 

26.22 

- 10.45 

<7 = 0.001 

1.000 

- 3 . 14 xl 0" 4 

2.756 

3 . 18 x 10 "* 

5.487 

- 4.27 

True 

1.000 

0 

2.760 

0 

5.400 

0 


After the model order is determined, the ERA algorithm estimates the modal frequencies 
and damping factors. Note that the true frequencies for this example are given in Eq. (21), and the 
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true damping factors are 0. The frequencies and damping factors estimated by the ERA method 
are given in Table (4), along with the true values. In constructing Table (4), we have chosen the 
number of modes as three in all cases, even though this is not clear from the singular values. The 
damping and frequency parameters in Table (4) clearly indicate that modes 2 and 3 have not been 
discerned from the noisier measurements. 

We now proceed to apply the combined ERA/MME algorithm to the same five sets of 
measurements used in the ERA algorithm. The assumed dynamic model varies from case to case. 
The results from Tables (3) and (4) were used to construct the models which were assumed for 
the MME algorithm. Thus, for the three highest noise levels, the MME used only the first mode 
identified by the ERA, and for the lower two noise levels, used the first two modes identified by the 
ERA. The MME algorithm produced state estimates for the measurement position. These esti- 
mates were then sampled at the original measurement times to produce “simulated” measurements 
which contain significantly less noise than the original measurements. Finally, the ERA algorithm 
is again applied, this time to the simulated measurements. Although a second application of this 
procedure may improve the reahzation/identification, as illustrated in Figure (1), we present re- 
sults for a single pass only. The singular values obtained by ERA processing of the sampled MME 
estimates are given in Table (5). These singular values indicate three modes for every noise level. 
The parameter identification results axe given in Table (6). All three frequencies are identified at 
all noise levels. The damping identification for the first two modes is also very good at all noise 
levels. 


Table 5. Singular values from the ERA/MME algorithm. 


a = 0.1 

<r = 0.05 

= 0.01 

<7 = 0.003 

a = 0.001 

II 

o 

o 

29.613 

29.496 

29.461 

29.460 

29.456 

29.456 

20.094 

20.048 

20.039 

20.030 

20.028 

20.028 

0.960 

1.091 

1.168 

1.2790 

1.264 

1.257 

0.706 

0.760 

0.788 

0.8610 

0.851 

0.845 

0.234 

0.175 

0.052 

0.0085 

0.0154 

0.026 

0.209 

0.160 

0.046 

0.0074 

0.0125 

0.022 

0.106 

0.078 

0.024 

0.0045 

0.0031 

10- 12 

0.075 

0.058 

0.021 

0.0031 

0.0025 

10- 12 

0.073 

0.057 

0.021 

0.0023 

0.0020 

icr 13 

0.072 

0.056 

0.020 

0.0014 

0.0014 

io-» 

0.069 

0.055 

0.020 

0.0013 

0.0013 

io - 13 

0.066 

0.050 

0.017 

0.0012 

0.0012 

10~ 13 


The results presented in Tables (5) and (6) indicate that the combined algorithm is capable 
of identifying modes with amplitudes as low as 1% of the noise. In each case, we have assumed 
the minimum model identified by the first pass of the ERA. Thus, for example, at a = 0.1, which 
is twice the amplitude of mode 2 and 100 times the amplitude of mode 3, the ERA algorithm 
identifies a single mode. Using only this one-mode model as input to the MME, the combined 
algorithm is still able to determine that the true model order is three, and give good accuracy in 
the parameter estimates. 
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Table 6. Parameter identification from the ERA/MME algorithm. 



Mode 1 

Mode 2 

Mode 3 

Noise 

Frequency 

Damping 

Frequency 

Damping 

Frequency 

Damping 

II 

o 

►“* 

0.996 

1.44xl0 -3 

— 

-4.09xl0“ 2 

5.408 

■■ 

<t =0.05 

0.998 

4.14x10“" 

I 

-7.97 xl0“ 3 

5.517 


a =0.01 

0.999 

-7.18xl0 _s 

2.768 

2.22xl0“ 3 

5.678 

-0.766 

<7 =0.003 

1.000 

-4.30xl0 -5 

2.7C6 

-1.62xl0“ 3 

5.589 

-0.142 

<7 =0.001 

1.000 

-4.61 xlO" 5 


2.00x10-" 

5.397 

0.044 

True 

1.000 

0 


0 

5.400 

0 


5.4 Sample Comparison With Extended Kalman Filter-Smoother 

To illustrate the potential advantages of the model error terms in the MME compared with 
process noise in filters, consider the following nonlinear problem. The truth is given by the equation 

x + t 2 

X = -- 

2x + t 

For illustration, the assumed model for the estimation algorithms is 

x = 0 (23) 

The measurements are perfect measurements of x . The MME is applied to this problem, and, 
for comparison, an extended Kalman filter- smoother is also used. The EKFS is modeled after 
the well-known Rauch- Tung- Streibel 21 filter-smoother, extended for the nonlinear problem. Since 
the assumed model is zero, the EKFS estimate must be constant between the measurements. The 
results of the two estimation approaches are shown in Figure (17). The model correction capability 
of the MME enables it to produce state estimates using a corrected version of the original model, 
so that the MME is not constant between measurements. The advantage of this approach is clear 
in Figure (17). Even though the measurements are perfect, the EKFS estimates between the 
measurements are poor. 



Summary and Conclusions 

In this paper, a new method for optimal post-experiment estimation has been described 
and its application demonstrated by numerous examples. The method is formulated to account 
for model error in a much more general and rigorous fashion than the process noise assumptions 
of typical filter algorithms. The state estimates are continuous and based on global measurement 
fits, compared with filter estimates which are discrete and based on local measurement fits. The 
MME method may give vastly improved state estimates when compared with filters for dynamic 
problems with significant model error, especially if the measurements are sparse and/or noisy. 

In the MME, model error is treated as an unknown and estimated along with the states. 
The estimated model error is automatically corrected in the original model in order to obtain the 
state estimates. For poorly modeled systems, this produces two significant benefits. First, the state 
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□ SMOOT MEO FILTER 
O PRESENT METHOO 

* measurements 


KALMAN FILTER TRUTH-M INUS-EST IMATE VARIANCE = 1.93 
PRESENT METHOO TRUTH-M INUS-EST IMATE VARIANCE = .052 
TUNEO • TRUTH-MEAS VAR = MEAS-EST VAR UO.-C-IOH 



Figure 17. Comparison of the MME with an extended Kalman filter-smoother in the 
absence of a model. The EKFS estimate must be constant between measurements, 
but the model correction in the MME shows clear advantages. 

estimates are based on a corrected model (unlike filters), and second, the model error estimates 
are available to aid in identification of an accurate model for subsequent use. 

Examples are given which demonstrate state estimation and exploitation of the model error 
estimates for both system identification and external force identification. The MME method shows 
considerable promise for use in numerous post-experiment estimation and identification problems, 
and should be considered for any application in which significant model error is suspected. 
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